Mesh Curvature and Geodesics¶

In this assignment, you will explore curvature on discrete surfaces. First you compute the Gaussian curvature and then the Mean curvature normal. Finally, you compute the shape operator for each vertex and use it to find both the mean curvature (again) and also the principal curvatures.

In [1]:
from pygel3d import hmesh, jupyter_display as jd
from numpy import zeros, dot
from numpy.linalg import norm
import numpy as np
from math import pi, acos, tan, acos, cos, sqrt
jd.set_export_mode(True)
In [2]:
# A couple of meshes are provided in the zip. Feel free to try more.
filename = "bunnygtest.obj"
# filename = "torus.obj"
is_bunny = ("bunny" in filename)
m = hmesh.obj_load(filename)
hmesh.triangulate(m)
jd.display(m)

Part 1: Compute and Integrate Gaussian Curvature¶

Compute the Gaussian curvature at all vertices. Visualize the Gaussian curvature on the mesh and compute the integral gaussian curvature over the mesh. Print out the result of the integral Gaussian curvature. Try this both for the bunny and the torus. Note that when visualizing the curvature, some vertices are likely to have extreme values. You probably need to clamp the values to a range in order to be able to get a meaningful visualization.

In [3]:
def clip_worst(a,fraction=0.1):
    a_sorted = np.sort(a)
    n = len(a)-1
    return np.clip(a, a_sorted[int(n*fraction)], a_sorted[int(n*(1.0-fraction))])


K = zeros((len(m.vertices()),))
pos = np.array(m.positions())
kint = 0
for v in m.vertices():
    # Insert code here which computes the Gaussian curvature for each vertex,
    # adds to the sum (kint), and assigns the value to K[v]
    
    # Your code here
    # ------------------

    phi_sum = 0
    area = 0
    
    for fid in m.circulate_vertex(v, mode="f"):
        # Sum all triangles areas
        area += m.area(fid)

        # Calculate phi with cosine rule
        vid_con = np.array(m.circulate_face(fid, mode="v"))

        v1,v2 = vid_con[vid_con != v]

        a = np.linalg.norm(pos[v1] - pos[v2])
        b = np.linalg.norm(pos[v] - pos[v1])
        c = np.linalg.norm(pos[v] - pos[v2])
        
        phi_sum += np.arccos((b*b + c*c - a*a)/(2*b*c))
    
    
    K[v] = (2*np.pi - phi_sum) / area  # Definition 8.3 
    
    kint += K[v] 
    
print("Integral Gaussian curvature", kint)

if is_bunny:
    K = clip_worst(K, fraction=0.1)

jd.display(m,data=K)
Integral Gaussian curvature 995864.8226074256

Part 2: Compute Mean Curvature Normal¶

In this exericse, you should compute the mean curvature normal $\mathbf{H} = \text{H}\cdot \mathbf{n}$ for each vertex in the mesh. On the mesh, we plot the length of this vector which is the mean curvature.

In [4]:
def normalize(v):
    return v / np.linalg.norm(v)

def vertex_areas(m):
    pos = m.positions()
    face_areas = [ m.area(f) for f in m.faces() ]
    ora = zeros(len(pos))
    for v in m.vertices():
        edges = [normalize(pos[vn]-pos[v]) for vn in m.circulate_vertex(v, mode='v')]
        areas = [face_areas[f] for f in m.circulate_vertex(v, mode='f')]
        indices = range(len(edges))
        angles = [acos(dot(edges[i],edges[(i+1)%len(edges)])) for i in indices]
        ora[v] = sum([areas[i]*angles[i]/pi for i in indices])
    return ora

def cot(theta):
    return np.cos(theta)/np.sin(theta)

def spatial_elem(m, v, hid):

    pos = m.positions()
    
    # alpha is the angle opposite to hid
    v1 = m.incident_vertex(hid)
    v2 = m.incident_vertex(m.next_halfedge(hid))

    a = np.linalg.norm(pos[v1] - pos[v2])
    b = np.linalg.norm(pos[v] - pos[v1])
    c = np.linalg.norm(pos[v] - pos[v2])

    alpha = np.arccos((a*a + c*c - b*b)/(2*a*c))

    # Beta is similar except for the opposite hid
    hid_opp = m.opposite_halfedge(hid)
    v2 = m.incident_vertex(m.next_halfedge(hid_opp))

    a = np.linalg.norm(pos[v1] - pos[v2])
    c = np.linalg.norm(pos[v] - pos[v2])

    beta = np.arccos((a*a + c*c - b*b)/(2*a*c))

    pi = pos[v]
    pj = pos[v1]

    return  (cot(alpha) + cot(beta)) * (pi - pj)


def mean_curvature_normal(m):
    N = m.no_allocated_vertices()
    Hn = np.zeros((N,3))
    
    # Insert code here which computes the Mean curvature normal H*n for each vertex 
    # using the contan formula.

    ora = vertex_areas(m)
    
    for v in m.vertices():
        
        # Your code here
        # ------------------

        spatial_sum = np.zeros(3)
        for hid in m.circulate_vertex(v, mode="h"):
            spatial_sum += spatial_elem(m, v, hid)
            
        Hn[v] = spatial_sum / (4.0 * ora[v])
        
    return Hn
In [5]:
# Display Mean Curvature
Hn = mean_curvature_normal(m)
if is_bunny:
    Hn = clip_worst(Hn, fraction = 0.2)
H = norm(Hn,axis=1)
print("Mean Curvature", H.shape)
jd.display(m, data=H)
Mean Curvature (3485,)

Part 3: Compute the Shape Operator¶

In this part, you compute the shape operator $\mathbf{S} = - \begin{bmatrix} a & b \\ b & c \end{bmatrix}$ and use it to find the principal curvatures $\text{kmin}$ and $\text{kmax}$ and the mean curvature $\text{H}$.

In [6]:
def shape_operator(m,v):
    pos = np.array(m.positions())
    
    # Insert code here which computes the shape operator matrix S that can be used to compute
    # the principal curvatures kmin and kmax and the mean curvature H
    
    # Your code here
    # ------------------


    # Create the a,b,n span
    n = np.array(m.vertex_normal(v))
    
    a_tilde_vec = np.random.rand(len(n)) # Create random array that is not parallel to the normal
    while np.abs(n @ a_tilde_vec.T) == 1:
        a_tilde_vec = np.random.rand(len(n))

    a_vec = (a_tilde_vec - n * (a_tilde_vec @ n))
    a_vec /= np.linalg.norm(a_vec)

    b_vec = np.cross(n,a_vec)

    # create U and F 
    v_neigh = list(m.circulate_vertex(v, mode="v"))
    p_neigh_centered_at_v = pos[v_neigh] - pos[v]

    T_T = np.vstack((a_vec, b_vec))
    uv = T_T @ p_neigh_centered_at_v.T
    U = np.vstack((uv[0]*uv[0]/2, uv[0]*uv[1], uv[1]*uv[1]/2)).T

    F = p_neigh_centered_at_v @ n

    a,b,c = (np.linalg.inv(U.T @ U) @ U.T) @ F

    # Create shape operator
    S = - np.array([[a,b],[b,c]])

    S_3d = T_T.T @ S @ T_T

    eig_vals = np.linalg.eigvals(S_3d)

    index_nearest_zero = np.argmin(np.abs(eig_vals))

    # Principal curvature
    kmin, kmax = sorted(eig_vals[np.arange(3) != index_nearest_zero]) 
    
    # Mean curvature
    H = (kmin + kmax) / 2

    return kmin, kmax, H
In [7]:
# Display Mean Curvature
H = zeros(len(m.vertices()))
for v in m.vertices():
    kmin, kmax, mean_curvature = shape_operator(m,v)
    H[v] = mean_curvature

H = clip_worst(H, fraction=0.1)
print("Mean Curvature")
jd.display(m,data=H)
Mean Curvature
In [8]:
# Display Principal Curvatures
Hn = mean_curvature_normal(m)
min_curvature = zeros(len(m.vertices()))
max_curvature = zeros(len(m.vertices()))
for v in m.vertices():
    kmin, kmax, mean_curvature = shape_operator(m,v)
    min_curvature[v] = kmin
    max_curvature[v] = kmax

min_curvature = clip_worst(min_curvature, fraction=0.1)
max_curvature = clip_worst(max_curvature, fraction=0.1)

print("Principal curvatures")
jd.display(m,data=min_curvature)
jd.display(m,data=max_curvature)
Principal curvatures

Questions¶

  • In the plot of Gaussian curvature, what differentiates regions with positive curvature from regions with negative curvature?

Answer: The positive curvatures are colored in green and the negative curvatures are colored in pink.

  • Would you expect the image to look different for mean curvature and, if so, how?

Answer: Yes, mean curvature would take the mean of the curvatre to all directions from a point. We would probably not see so obviously the sharp changes between the inner and the outer parts.

  • Explain the value you got for the integral Gaussian curvature. Why is it different for the bunny and the torus?

Answer: The value gives us information about the shape of the object. This is basically the sum of the total curvature of the object. The bunny is all just curving "out" so the curvatures are positive in most points so the value is big and positive. The Torus has hole and many points have negative curvatures (also many points with positive). So the value for the torus is much closer to zero.

  • On the torus, why is the Gaussian curvature positive and negative, but the mean curvature always positive?

Answer: The outside of the Torus is a positive curvature and the inside of the torus is a negative curvature. The mean curvature is always positive is the average of both the positive and the negative, but the outside part influences the mean more and therefore the mean is always positive for the torus.

  • On the torus, the maximum principal curvature is almost constant, but the minimum curvature varies a lot. Why is that?

Answer: On a torus, the maximum curvature stays the same because it's like the thickness of the torus, that doesn't change. The minimum curvature changes a lot because it is like going from the inside of the hole to the outside, and that shape changes a lot.